⚡ Pluto.jl ⚡

xxxxxxxxxx
 
begin
    using Plots
    using PlutoUI
    using Statistics
end

8.4 s

Deriving a SIR epidemiological model

In this notebook we will go from a discrete epidemiological model to its continuous counterpart. We will go from a (discrete) system of difference equiations to a system of (continuous) differential equations. But first: SIR.

Susceptible, Infected, Recovered (SIR)

In this workhorse of epidemiology we model agents as belonging to a certain class, or compartment. Either they are

susceptible to get infected,
infected, or
removed/recovered

The models you may have heard about during the COVID19 pandemic are more sophisticated versions of what we will see here, but they share some important common features. At the end of this notebook, we will want to generate the following plot. It shows the number of agents belonging to each class at a certain point in time after the outbreak of the disease:

4.0 ms

xxxxxxxxxx
 
discrete_time_SIR_plot

37.0 ns

We will make this very simple. Instead of trying to model in great detail which type of agent (old, young, where, what time of the day, doing which activity etc) is likely to get infected, we will try to condense the overall probability of infection into a single number. Let's start to model the recovery from an infection $I \to R$ :

we have $N$ people infected at time $0$ . If each has probability $p$ to recover each day, how many are still infected at day number $t$ ?

For any given person in $I$ state, the probability to recover in a certain time step (i.e. each day) is a

3.2 μs

bernoulli (generic function with 1 method)

xxxxxxxxxx
 
bernoulli(p) = rand() < p 

21.6 μs

xxxxxxxxxx
 
N = 200

27.0 ns

step! (generic function with 1 method)

xxxxxxxxxx
 
function step!(infectious, p)
    for i in 1:length(infectious)
        
        if infectious[i] && bernoulli(p)
            infectious[i] = false
        end
    end
    
    return infectious
end

28.1 μs

xxxxxxxxxx
 
T = 100

40.0 ns

BitArray{1}1BitArray{1}: 1

true

true

true

true

true

true

true

true

true

true

true

true

true

true

true

true

true

true

true

true

more191

true

192

true

193

true

194

true

195

true

196

true

197

true

198

true

199

true

200

true

2BitArray{1}: 1

true

true

true

true

true

true

true

true

true

true

true

true

true

true

true

true

true

true

true

true

more191

true

192

true

193

true

194

true

195

true

196

true

197

true

198

true

199

true

200

true

3BitArray{1}: 1

true

true

true

true

true

false

true

true

true

true

false

true

true

true

false

true

true

true

true

true

more191

true

192

true

193

true

194

true

195

true

196

true

197

true

198

true

199

true

200

true

4BitArray{1}: 1

true

true

true

true

true

false

false

false

true

true

false

true

true

true

false

true

true

true

true

true

more191

true

192

true

193

true

194

true

195

true

196

true

197

true

198

true

199

true

200

true

5BitArray{1}: 1

true

true

true

true

true

false

false

false

true

true

false

false

true

true

false

true

true

true

true

true

more191

true

192

true

193

true

194

true

195

true

196

true

197

true

198

true

199

true

200

true

6BitArray{1}: 1

true

true

true

true

true

false

false

false

true

true

false

false

true

true

false

true

true

true

true

true

more191

true

192

false

193

true

194

true

195

true

196

true

197

true

198

true

199

true

200

true

7BitArray{1}: 1

true

true

true

true

true

false

false

false

true

true

false

false

true

true

false

true

true

true

true

true

more191

true

192

false

193

true

194

true

195

false

196

true

197

true

198

true

199

true

200

true

8BitArray{1}: 1

true

true

true

true

true

false

false

false

true

true

false

false

true

true

false

true

true

true

true

true

more191

true

192

false

193

true

194

true

195

false

196

true

197

true

198

false

199

true

200

true

9BitArray{1}: 1

true

true

true

true

true

false

false

false

true

true

false

false

true

true

false

true

true

true

true

true

more191

true

192

false

193

true

194

true

195

false

196

true

197

false

198

false

199

true

200

true

10BitArray{1}: 1

true

false

true

true

true

false

false

false

true

true

false

false

false

true

false

true

true

true

true

true

more191

true

192

false

193

true

194

true

195

false

196

true

197

false

198

false

199

true

200

true

11BitArray{1}: 1

true

false

true

true

true

false

false

false

true

true

false

false

false

true

false

true

true

true

true

true

more191

true

192

false

193

true

194

true

195

false

196

true

197

false

198

false

199

false

200

true

12BitArray{1}: 1

true

false

true

true

true

false

false

false

true

true

false

false

false

true

false

true

false

true

true

true

more191

true

192

false

193

false

194

true

195

false

196

true

197

false

198

false

199

false

200

true

13BitArray{1}: 1

true

false

true

true

true

false

false

false

true

true

false

false

false

true

false

true

false

true

true

true

more191

true

192

false

193

false

194

true

195

false

196

true

197

false

198

false

199

false

200

false

14BitArray{1}: 1

true

false

true

false

true

false

false

false

true

true

false

false

false

true

false

true

false

true

true

true

more191

true

192

false

193

false

194

true

195

false

196

true

197

false

198

false

199

false

200

false

15BitArray{1}: 1

true

false

true

false

true

false

false

false

true

true

false

false

false

true

false

true

false

true

true

true

more191

true

192

false

193

false

194

true

195

false

196

true

197

false

198

false

199

false

200

false

16BitArray{1}: 1

true

false

true

false

true

false

false

false

true

true

false

false

false

true

false

true

false

true

true

true

more191

true

192

false

193

false

194

true

195

false

196

true

197

false

198

false

199

false

200

false

17BitArray{1}: 1

true

false

true

false

true

false

false

false

true

true

false

false

false

true

false

true

false

true

true

true

more191

true

192

false

193

false

194

true

195

false

196

true

197

false

198

false

199

false

200

false

18BitArray{1}: 1

true

false

true

false

true

false

false

false

true

true

false

false

false

true

false

true

false

true

true

true

more191

true

192

false

193

false

194

true

195

false

196

true

197

false

198

false

199

false

200

false

19BitArray{1}: 1

true

false

true

false

true

false

false

false

true

true

false

false

false

false

false

true

false

true

true

true

more191

true

192

false

193

false

194

true

195

false

196

true

197

false

198

false

199

false

200

false

20BitArray{1}: 1

true

false

true

false

true

false

false

false

true

true

false

false

false

false

false

true

false

false

true

true

more191

true

192

false

193

false

194

true

195

false

196

true

197

false

198

false

199

false

200

false

more92BitArray{1}: 1

false

false

false

false

false

false

false

false

false

false

false

false

false

false

false

false

false

false

false

false

more191

false

192

false

193

false

194

true

195

false

196

false

197

false

198

false

199

false

200

false

93BitArray{1}: 1

false

false

false

false

false

false

false

false

false

false

false

false

false

false

false

false

false

false

false

false

more191

false

192

false

193

false

194

true

195

false

196

false

197

false

198

false

199

false

200

false

94BitArray{1}: 1

false

false

false

false

false

false

false

false

false

false

false

false

false

false

false

false

false

false

false

false

more191

false

192

false

193

false

194

true

195

false

196

false

197

false

198

false

199

false

200

false

95BitArray{1}: 1

false

false

false

false

false

false

false

false

false

false

false

false

false

false

false

false

false

false

false

false

more191

false

192

false

193

false

194

true

195

false

196

false

197

false

198

false

199

false

200

false

96BitArray{1}: 1

false

false

false

false

false

false

false

false

false

false

false

false

false

false

false

false

false

false

false

false

more191

false

192

false

193

false

194

false

195

false

196

false

197

false

198

false

199

false

200

false

97BitArray{1}: 1

false

false

false

false

false

false

false

false

false

false

false

false

false

false

false

false

false

false

false

false

more191

false

192

false

193

false

194

false

195

false

196

false

197

false

198

false

199

false

200

false

98BitArray{1}: 1

false

false

false

false

false

false

false

false

false

false

false

false

false

false

false

false

false

false

false

false

more191

false

192

false

193

false

194

false

195

false

196

false

197

false

198

false

199

false

200

false

99BitArray{1}: 1

false

false

false

false

false

false

false

false

false

false

false

false

false

false

false

false

false

false

false

false

more191

false

192

false

193

false

194

false

195

false

196

false

197

false

198

false

199

false

200

false

100BitArray{1}: 1

false

false

false

false

false

false

false

false

false

false

false

false

false

false

false

false

false

false

false

false

more191

false

192

false

193

false

194

false

195

false

196

false

197

false

198

false

199

false

200

false

101BitArray{1}: 1

false

false

false

false

false

false

false

false

false

false

false

false

false

false

false

false

false

false

false

false

more191

false

192

false

193

false

194

false

195

false

196

false

197

false

198

false

199

false

200

false

xxxxxxxxxx
 
# try this out: all are infected on day one.
# let's record how many are in each state I or R after i days
begin
    infected = trues(N)  # all infected
        
    results = [copy(step!(infected, 0.05)) for i in 1:T]
    # we `copy` here because we keep modifying the same array throughout
    pushfirst!(results, trues(N))
end

83.3 μs

xxxxxxxxxx
 
@bind i Slider(1:T, show_value=true)

23.8 ms

5.8 s

ok great, now we know how many are $I$ or $R$ after $t$ days with a certain probability $p$ . Let's just condense this plot above and count how many infected there are each day!

8.3 μs

simulate_recovery (generic function with 1 method)

xxxxxxxxxx
 
function simulate_recovery(p, T)
    infectious = trues(N)
    num_infectious = [N]
    
    for t in 1:T
        step!(infectious, p)
        push!(num_infectious, count(infectious))
    end
    
    return num_infectious
end

31.7 μs

xxxxxxxxxx
 
begin
    pp = 0.05
    
    plot(simulate_recovery(pp, T), label="run 1", alpha=0.5, lw=2, m=:o)
    plot!(simulate_recovery(pp, T), label="run 2", alpha=0.5, lw=2, m=:o)
    
    xlabel!("time t")
    ylabel!("number infectious")
end

652 ms

xxxxxxxxxx
 
all_data = [simulate_recovery(pp, T) for i in 1:30];  # do 30 runs

25.6 ms

515 ms

279 ms

xxxxxxxxxx
 
begin
    plot(replace.(all_data, 0.0 => NaN), 
        yscale=:log10, alpha=0.3, leg=false, m=:o, ms=1,
        size=(500, 400))
    
    plot!(mean(all_data), yscale=:log10, lw=3, c=:red, m=:o, label="mean", alpha=0.5)
    
    xlabel!("time t")
    ylabel!("number infectious")
end
​
​

372 ms

Deterministic dynamics for the mean: Intuitive derivation

20.5 μs

The mean seems to behave in a rather predictable way over time. Can we derive this?

Let $I_{t}$ be the number of infectious people at time $t$ . This decreases because some people recover. Since people recover with probability $p$ , the number of people that recover at time $t$ is, on average, $p I_{t}$ . [Note that one time unit corresponds to one sweep of the simulation.]

$I_{t + 1} = I_{t} - p I_{t}$

$I_{t + 1} = I_{t} (1 - p) .$

8.5 μs

At time $t$ there are $I_{t}$ infectious. How many decay? Each decays with probability $p$ , so on average $p I_{t}$ recover, so are removed from the number of infectious, giving the change

$Δ I_{t} = I_{t + 1} - I_{t} = - p I_{t}$

8.9 μs

We can rearrange and solve the recurrence relation:

$I_{t + 1} = (1 - p) I_{t} .$ so $I_{t + 1} = (1 - p) (1 - p) I_{t - 1} = (1 - p)^{2} I_{t - 1}$

and hence solve the recurrence relation:

$I_{t} = (1 - p)^{t} I_{0} .$

16.8 μs

Let's compare the exact and numerical results:

8.6 μs

exact

Float641

200.0

190.0

180.5

171.475

162.901

154.756

147.018

139.667

132.684

126.05

119.747

113.76

108.072

102.668

97.535

92.6582

88.0253

83.6241

79.4429

75.4707

more92

1.87879

1.78485

1.69561

1.61083

1.53029

1.45377

1.38108

1.31203

100

1.24643

101

1.18411

5.4 μs

xxxxxxxxxx
 
begin
    plot(mean(all_data), m=:o, label="numerical mean")
    plot!(exact, lw=3, label="analytical mean")
end
    

65.1 ms

They agree well, as they should. The agreement is expected to be better (i.e. the fluctuations smaller) for a larger population.

8.7 μs

Continuous time

If we look at the graph of the mean as a function of time, it seems to follow a smooth curve. Indeed it makes sense to ask not only how many people have recovered each day, but to aim for finer granularity.

Suppose we instead increment time in steps of $δ t$ ; the above analysis was for $δ t = 1$ .

Then we will need to adjust the probability of recovery in each time step. It turns out that to make sense in the limit $δ t \to 0$ , we need to choose the probability $p (δ t)$ to recover in time $t$ to be proportional to $δ t$ :

$p (δ t) ≃ λ δ t,$

where $λ$ is the recovery rate. Note that a rate is a probability per unit time.

We get

17.2 μs

$I (t + δ t) - I (t) ≃ - λ δ t I (t)$

2.5 μs

Dividing by $δ t$ gives

$\frac{I (t + δ t) - I (t)}{δ t} ≃ - λ I (t)$

We recognise the left-hand side as the definition of the derivative when $δ t \to 0$ . Taking that limit finally gives

9.7 μs

$\frac{d I (t)}{d t} = - λ I (t)$

That is, we obtain an ordinary differential equation that gives the solution implicitly. We know this has a general solution of the form

$I (t) = C \exp (- λ t),$

and together with an initial value of $I (0) = I_{0}$ we get (just set $t = 0$ and use the given value)

11.6 μs

$I (t) = I_{0} \exp (- λ t) .$

Starting at time zero ( $t = 0$ ) and at a number of $I_{0}$ infected, this is exponential decay at rate $λ$ . Of course with a positive number (i.e. no minus sign, as here), we get exponential growth.

9.8 μs

SIR model

3.1 μs

Now let's extend the procedure to the full SIR model, $S \to I \to R$ . Since we already know how to deal with recovery, consider just the SI model, where susceptible agents are infected via contact, with a certain probability.

xxxxxxxxxx
 
md"""
Now let's extend the procedure to the full SIR model, $S \to I \to R$. Since we already know how to deal with recovery, consider just the SI model, where susceptible agents are infected via contact, with a certain probability.
"""

3.1 μs

Let's denote by $S_{t}$ and $I_{t}$ be the number of susceptible and infectious people at time $t$ , respectively, and by $N$ the total number of people.

On average, in each sweep each infectious individual has the chance to interact with one other individual. That individual is chosen uniformly at random from the total population of size $N$ . But a new infection occurs only if that chosen individual is susceptible, which happens with probability $S_{t} / N$ . Then, upon meeting a susceptible person, the infection is transmitted only with a certain probability, $b$ say.

Hence the change in the number of infectious people after that step is:

15.2 μs

$Δ I_{t} = I_{t + 1} - I_{t} = b I_{t} (\frac{S_{t}}{N})$

3.1 μs

It is useful to normalize by $N$ , so we define

$s_{t} := \frac{S_{t}}{N}; i_{t} := \frac{I_{t}}{N}; r_{t} := \frac{R_{t}}{N}$

4.3 μs

Including recovery with probability $c$ we obtain the discrete-time SIR model:

5.1 μs

$\begin{aligned} s_{t + 1} & = s_{t} - b s_{t} i_{t} \\ i_{t + 1} & = i_{t} + b s_{t} i_{t} - c i_{t} \\ r_{t + 1} & = r_{t} + c i_{t} \end{aligned}$

xxxxxxxxxx
 
md"""
$$\begin{align}
s_{t+1} &= s_t - b \, s_t \, i_t \\
i_{t+1} &= i_t + b \, s_t \, i_t - c \, i_t\\
r_{t+1} &= r_t + c \, i_t
\end{align}$$
"""

1.6 μs

And again we can allow the processes to occur in steps of length $δ t$ and take the limit $δ t \to 0$ . With rates $β$ and $γ$ we obtain the standard (continuous-time) SIR model:

3.6 μs

$\begin{aligned} \frac{d s (t)}{d t} & = - β s (t) i (t) \\ \frac{d i (t)}{d t} & = + β s (t) i (t) & - γ i (t) \\ \frac{d r (t)}{d t} & = & + γ i (t) \end{aligned}$

xxxxxxxxxx
 
md"""
$$\begin{align}
\textstyle \frac{ds(t)}{dt} &= -\beta \, s(t) \, i(t) \\
\textstyle \frac{di(t)}{dt} &= +\beta \, s(t) \, i(t) &- \gamma \, i(t)\\
\textstyle \frac{dr(t)}{dt} &= &+ \gamma \, i(t)
\end{align}$$
"""

1.7 μs

Note that no analytical solutions of these (simple) nonlinear ODEs are known as a function of time!

2.5 μs

Below is a simulation of the discrete-time model. Note that the simplest numerical method to solve (approximately) the system of ODEs, the Euler method, basically reduces to solving the discrete-time model! A whole suite of more advanced ODE solvers is provided in the Julia DiffEq ecosystem. We will in another notebook introduce that package, together with a continuous version of the SIR model.

7.8 μs

xxxxxxxxxx
 
begin
    ts = 1:length(SIR)
    discrete_time_SIR_plot = plot(ts, [x.s for x in SIR], 
        m=:o, label="S", alpha=0.2, linecolor=:blue, leg=:right, size=(400, 300))
    plot!(ts, [x.i for x in SIR], m=:o, label="I", alpha=0.2)
    plot!(ts, [x.r for x in SIR], m=:o, label="R", alpha=0.2)
    
    xlims!(0, 500)
end

579 ms

xxxxxxxxxx
 
begin
    NN = 100
    
    SS = NN - 1
    II = 1
    RR = 0
end

38.0 ns

0.99

0.01

0.0

xxxxxxxxxx
 
ss, ii, rr = SS/NN, II/NN, RR/NN

31.0 ns

0.1

0.01

xxxxxxxxxx
 
p_infection, p_recovery = 0.1, 0.01

34.0 ns

discrete_SIR (generic function with 2 methods)

xxxxxxxxxx
 
function discrete_SIR(s0, i0, r0, T=1000)
​
    s, i, r = s0, i0, r0
    
    results = [(s=s, i=i, r=r)]
    
    for t in 1:T
​
        Δi = p_infection * s * i
        Δr = p_recovery * i
        
        s_new = s - Δi
        i_new = i + Δi - Δr
        r_new = r      + Δr
​
        push!(results, (s=s_new, i=i_new, r=r_new))
​
        s, i, r = s_new, i_new, r_new
    end
    
    return results
end

46.3 μs

SIR

NamedTuple{(:s, :i, :r),Tuple{Float64,Float64,Float64}}1s

0.99

0.01

0.0

0.98901

0.01089

0.0001

0.987933

0.0118581

0.0002089

0.986761

0.0129111

0.000327481

0.985487

0.014056

0.000456592

0.984102

0.0153006

0.000597151

0.982597

0.0166533

0.000750157

0.98096

0.0181231

0.000916691

0.979182

0.0197197

0.00109792

10s

0.977251

0.0214534

0.00129512

11s

0.975155

0.0233354

0.00150965

12s

0.972879

0.0253777

0.00174301

13s

0.97041

0.0275928

0.00199678

14s

0.967733

0.0299945

0.00227271

15s

0.96483

0.0325973

0.00257266

16s

0.961685

0.0354164

0.00289863

17s

0.958279

0.0384681

0.00325279

18s

0.954593

0.0417698

0.00363748

19s

0.950605

0.0453394

0.00405517

20s

0.946295

0.049196

0.00450857

more992s

3.65235e-5

8.28801e-5

0.999881

993s

3.65232e-5

8.20516e-5

0.999881

994s

3.65229e-5

8.12314e-5

0.999882

995s

3.65226e-5

8.04194e-5

0.999883

996s

3.65223e-5

7.96155e-5

0.999884

997s

3.65221e-5

7.88196e-5

0.999885

998s

3.65218e-5

7.80317e-5

0.999885

999s

3.65215e-5

7.72517e-5

0.999886

1000s

3.65212e-5

7.64794e-5

0.999887

1001s

3.65209e-5

7.57149e-5

0.999888

xxxxxxxxxx
 
SIR = discrete_SIR(ss, ii, rr)

11.5 ms